1 Descripción del ejercicio

El presente ejercicio se engloba dentro de la materia de Sistemas de Información Geográfica y Ecología Espacial: Aplicaciones del Máster Geoforest de la Universidad de Córdoba.

El trabajo consiste en aplicar las herramientas aprendidas en clase, con lo que los productos finales esperados son:

1.1 Datos de partida

Para la realización del ejercicio partimos de los siguientes datos:

  • limite_yunquera.shp: límite de la zona de estudio (EPSG: 25830).

  • vegetacion_yunquera.shp: mapa de vegetación de la zona de estudio.

  • inventario_pinsapo.shp: puntos de inventario de parcelas de pinsapo (Abies pinsapo).

  • PNOA_MDT05_ETRS89_HU30_1051_LID.tif: modelo digital del terreno con tamaño de píxel de 5x5 metros.

Estos datos se encuentran en la carpeta Cartografia/CartografiaOriginal.

El resto de datos de partida se obtendrán a lo largo de la resolución del ejercicio.

1.2 Zona de estudio

Antes de comenzar a resolver el ejercicio conviene realizar una breve descripción de la zona de estudio.

Yunquera es un municipio perteneciente a la provincia de Málaga (Andalucía), que se sitúa al oeste de la provincia.

Al oeste del municipio de Yunquera, se encuentra el monte Pinar de Yunquera. En este monte encontramos especies arbóreas como Abies pinsapo, Castanea sativa, Cedrus atlantica, Quercus faginea, Quercus ilex subsp. ballota, Pinus halepensis, Pinus nigra y Olea europaea var. europaea. Estos se encuentran formando distintas formaciones entre las que destacan los pinsapares en masas monoespecíficas y pluriespecíficas.

En el siguiente webmap se pueden explorar las diferentes formaciones vegetales presentes en el área de estudio (hacer click en las parcelas para explorar).

yunquera <- read_sf(here('Cartografia/CartografiaOriginal/limite_yunquera.shp'),
                    crs = 25830)
vegetacion <- read_sf(here('Cartografia/CartografiaOriginal/vegetacion_yunquera.shp'),
                    crs = 25830)


2 Resolución del ejercicio

Se ha dividido el ejercicio en 5 bloques de acuerdo a los objetivos del mismo.

2.1 Flujo de trabajo

Se presenta el flujo de trabajo como iframe para poder hacer zoom o abrir el diagrama en diagrams.net. También se puede visualizar directamente aquí (recomendado).


2.2 Mapa de biodiversidad

En esta sección se construirá un mapa donde se representará el Índice de Shanon:

\[ \tag1 H' = - \sum^s_{i=1} p_i \times log_2 p_i \]


donde:

  • i: cada una de las distintas especies

  • s: total de especies distintas

  • \(p_i\): abundancia relativa de cada especie en la comunidad (N individuos de la especie i / N total de individuos).

2.2.1 Descarga de ocurrencias

En primer lugar debemos descargar las ocurrencias de GBIF. En este caso se ha creado un rectángulo en la plataforma web y se han extraído todas las ocurrencias que se encontraban dentro de este rectángulo (correspondiente a la zona de estudio).

A continuación se importan las ocurrencias en R con el nombre de puntos. El siguiente paso es convertir estos puntos en un objeto espacial. Para ello se utiliza la función sf::st_as_sf(). La salida de esta operación (puntos_sf) se transforma al sistema de coordenadas en el que se está trabajando (EPSG: 25830) llamando al objeto ocurrencias. En la Fig. 1 comprobamos que hemos elegido correctamente el rectángulo, y que nuestros puntos se extienden más allá de la zona de estudio asegurando la inclusión de todas las ocurrencias.

puntos <- read.csv(file = here('Cartografia/Tablas/ocurrences.csv'),
                   sep = "\t",
                   quote = "") |> 
  dplyr::select(scientificName,species,decimalLatitude,decimalLongitude) |> 
  drop_na()

puntos_sf <- st_as_sf(puntos,
                      coords = c('decimalLongitude','decimalLatitude'),
                      crs = 4326)

puntos_sf <- st_transform(puntos_sf,
                            crs = 25830)
<center>Fig. 1. Localización de los puntos descargados de GBIF con respecto al área de estudio</center>

Fig. 1. Localización de los puntos descargados de GBIF con respecto al área de estudio


El siguiente paso será delimitar los puntos escogidos al rectánculo mínimo envolvente (bounding box) del área de estudio mediante la función sf::st_crop(). La razón de escoger todos los puntos incluidos dentro del rectángulo y no solamente los que están dentro del área de estudio es que como veremos en la siguiente sección trabajaremos con comunidades ecológicas de forma cuadrangular.

ocurrencias <- st_crop(puntos_sf, yunquera)

Podemos ver el resultado en la Fig. 2, que ascienden a un total de 16906 ocurrencias correspondientes a 1465 especies diferentes.

<center>Fig. 2. Ocurrencias de GBIF dentro del monte de Yunquera</center>

Fig. 2. Ocurrencias de GBIF dentro del monte de Yunquera


2.2.2 Creación de comunidades ecológicas

Las comunidades ecológicas se asumirán que son los polígonos de la capa de vegetaciones. Para ello, en primer lugar, se crea un identificador único llamado ID y se seleccionan las columnas que más nos interesan. En el siguiente paso se crea una unión (st_join()) entre las ocurrencias y el mapa de vegetación según su intersección. De este modo podemos asignar el ID de la comunidad a cada ocurrencia. Finalmente se transforman las columnas necesarias en data.frame para los posteriores cálculos.

# Crear un ID único
vegetacion <- vegetacion |> 
  cbind(data.frame(ID = sprintf(paste("GID%0",
                                      nchar(nrow(vegetacion)),
                                      "d",
                                      sep=""), 
                                1:nrow(vegetacion)))) |> 
  dplyr::select(ID, D_ARBO_PRE:D_SUEL_PRE, D_USO:COMENTARIO)

# Añadir ID de comunidad ecologica a cada ocurrencia
ocurrencias_id <- ocurrencias |> 
  st_join(vegetacion, join = st_intersects) |> 
  mutate(id_comunidad = factor(ID),
         scientific = scientificName) |> 
  dplyr::select(-ID,-scientificName)

# Crear data frame 
bio <- data.frame(id_comunidad = ocurrencias_id$id_comunidad,
                  scientific = ocurrencias_id$scientific) |> 
  drop_na()
Tabla 1. Encabezado de la tabla de datos de ocurrencias en las distintas comunidades
id_comunidad scientific
GID024 Phoenicurus ochruros (S.G.Gmelin, 1774)
GID031 Abies pinsapo Boiss.
GID043 Erinacea anthyllis Link
GID018 Erinacea anthyllis Link
GID075 Ophrys lutea Cav.
GID100 Abies pinsapo Boiss.


En la Fig. 3 tenemos una representación gráfica de como quedarían estas comunidades con respecto a las ocurrencias, y en la Fig. 4 vemos la densidad de puntos en la zona.

<center>Fig. 3. Comunidades ecológicas en el monte de Yunquera. Los puntos son todas las ocurrencias consideradas en los análisis posteriores</center>

Fig. 3. Comunidades ecológicas en el monte de Yunquera. Los puntos son todas las ocurrencias consideradas en los análisis posteriores


<center>Fig. 4. Mapa de densidad de ocurrencias en el monte de Yunquera</center>

Fig. 4. Mapa de densidad de ocurrencias en el monte de Yunquera


2.2.3 Cálculo del índice de Shannon

El último paso de esta sección será calcular el Índice de Shannon (H). Esto lo haremos en cuatro pasos:

  1. Calcular el número de especies por comunidad. En el siguiente cuadro de código se agrupa por comunidad y nombre (group_by), se genera una columna con el número de especies (count()), y se le cambia el nombre a la columna. En la Tabla 2 vemos el encabezado del resultado.
T_num_ind_sp_com <- bio |> 
  group_by(id_comunidad, scientific) |> 
  count() |> 
  rename(num_ind_sp_com = n)
Tabla 2. Encabezado de la tabla del número de especies por comunidad ecológica
id_comunidad scientific num_ind_sp_com
GID001 Avena sativa L. 2
GID001 Thymus mastichina subsp. mastichina 9
GID001 Ulex baeticus subsp. baeticus 5
GID002 Abies pinsapo Boiss. 14
GID002 Asplenium trichomanes L. 1
GID002 Juniperus phoenicea L. 2


  1. Calcular el número de ocurrencias por comunidad. Realizamos un procedimiento similar, pero en este caso solamente agrupamos por comunidad. Vemos el encabezado en la Tabla 3.
T_num_ind_com <- bio |> 
  group_by(id_comunidad) |> 
  count() |> 
  rename(num_ind_com = n)
Tabla 3. Encabezado de la tabla del número de ocurrencias por comunidad ecológica
id_comunidad num_ind_com
GID001 16
GID002 40
GID004 42
GID005 26
GID006 2
GID007 29


  1. El siguiente paso calcula el Índice de Shannon (T_Shannon). Para ello se unen las tablas haciendo un left_join() que incluye todas las columnas de la Tabla 1 de forma que cada comunidad tendrá asociada sus especies, el número de ocurrencias de cada una, y el número total de ocurrencias de la cuadrícula de esa especie. A continuación se crean los coeficientes \(p_i\) y \(\log p_i\) que se utilizarán para el cálculo del índice. En la Tabla 4 vemos la estructura de los datos hasta este punto.
    Finalmente, se calcula el índice de Shannon, se seleccionan solamente las columnas necesarias y se eliminan los duplicados.
T_Shannon <- T_num_ind_sp_com |> 
  left_join(T_num_ind_com,
            by = 'id_comunidad') |> 
  mutate(pi = num_ind_sp_com/num_ind_com,
         logpi = log2(pi)) |> 
  group_by(id_comunidad) |> 
  mutate(shannon = sum(logpi * pi)*-1) |> 
  dplyr::select(id_comunidad, shannon) |> 
  distinct()
Tabla 4. Coeficientes para el cálculo del Índice de Shannon
id_comunidad scientific num_ind_sp_com num_ind_com pi logpi
GID001 Avena sativa L. 2 16 0.1250 -3.000000
GID001 Thymus mastichina subsp. mastichina 9 16 0.5625 -0.830075
GID001 Ulex baeticus subsp. baeticus 5 16 0.3125 -1.678072
GID002 Abies pinsapo Boiss. 14 40 0.3500 -1.514573
GID002 Asplenium trichomanes L. 1 40 0.0250 -5.321928
GID002 Juniperus phoenicea L. 2 40 0.0500 -4.321928


  1. El último paso consiste en dar el valor de Índice de Shannon a cada una de las comunidades. Se añade un paso de reemplazar por 0 en caso de existir NA. Vemos en el webmap la representación final.
# Valor H a cada cuadricula
SF_shannon <- vegetacion |> 
  left_join(y = T_Shannon,
            by = c("ID" = "id_comunidad"))

# Sustituir NA por 0
SF_shannon$shannon <- replace_na(SF_shannon$shannon, 0)


2.3 Segmentación

La segmentación consiste en dividir el territorio en unidades homogéneas en cuanto a unas determinadas características de las variables de entrada. El procedimiento consiste de los siguientes pasos:

  1. Elegir las variables de entrada

  2. Convertir a formato raster y homogeneizar las escalas

  3. Construir un raster virtual

  4. Segmentación del territorio

2.3.1 Variables de entrada

En esta sección se realizarán los pasos 1 y 2 previos. Para ello, construiremos un total de 4 variables (ver cuadro siguiente).

Variables de entrada

Densidad de pinsapo: información sobre la vegetación

Orientaciones: información sobre el estrés abiótico derivado de la situación de solana y umbría

Distancia a canales: información sobre zonas de posible compensación hídrica

Índice de Shannon: información sobre la biodiversidad


Densidad de pinsapo

En primer lugar cargamos la capa, la cual se encuentra en CRS 25830, y la visualizamos con los límites del monte para comprobar que todo es correcto.

SF_pinsapo <- read_sf(here("Cartografia/CartografiaOriginal/inventario_pinsapo.shp"),
                      crs = 25830)


El siguiente paso será crear un raster de densidad. En ese sentido, se utilizará la interpolación inversa a la distancia (IDW).

  1. Primero se transforman los puntos a WGS 1984 ya que las funciones siguientes trabajan este sistema de coordenadas.
SF_pinsapo84 <- st_transform(SF_pinsapo, crs = 4326)
  1. A continuación obtenemos el rectángulo mínimo envolvente(BBox84). Se da una toleracia de 0.005º dado que al reproyectar la imagen se pierde información.
BBox84 <- st_bbox(SF_pinsapo84)


BBox84[c(1,3)] <- BBox84[c(1,3)] + c(-0.005,0)
BBox84[c(2,4)] <- BBox84[c(2,4)] + c(-0.005,0)
  1. Creamos un patrón de puntos con la función spatstat::ppp() utilizando el vector de coordendas X e Y, marks es la variable que utilizaremos para interpolar (Nparc), y en un objeto owin() ponemos las coordendas X e Y del rectángulo mínimo envolvente, que será la extensión de salida.
ppp_pinsapo <- ppp(st_coordinates(SF_pinsapo84)[,1],
                   st_coordinates(SF_pinsapo84)[,2],
                   marks = SF_pinsapo$Nparc,
                   window = owin(BBox84[c(1,3)],
                                 BBox84[c(2,4)]))
  1. Con la función spatstat::idw() realizamos la interpolación y se transforma a raster, que en este caso a diferencia de QGIS no tenemos opción de escoger el número de vecinos de cuáles hacer la interpolación (no obstante, la diferencia del raster de salida es mínima).
R_pinsapo <- ppp_pinsapo |> 
  idw(power = 2, at = 'pixels') |> 
  raster()
  1. Se asigna la proyección y se reproyecta al sistema 25830 cambiando la resolución a 10 metros.
# Asignar proyección
crs(R_pinsapo) <- crs(SF_pinsapo84)

# Transformar a 25830 y asignar resolución espacial
R_pinsapo <- projectRaster(R_pinsapo,
                           crs = 25830,
                           res = 10)
  1. Finalmente se hace una reclasificación de los valores del ráster y se recorta al área de estudio.
# Matriz de entrada para reclasificación
classMatrix <- matrix(c(0,50,0,
                        50,200,100,
                        200,400,300,
                        400,600,500,
                        600,800,700,
                        800,Inf,900),
                      byrow = T,
                      ncol = 3)

# Reclasificación y recorte al área de estudio
R_pinsapo <- R_pinsapo |> 
  reclassify(rcl = classMatrix) |> 
  mask(yunquera)

En el siguiente webmap podemos ver el resultado de los análisis anteriores.


Orientaciones

La segunda de las capas se corresponde con la de orientaciones. En este caso cargamos el modelo digital del terreno (mdt), calculamos las orientaciones en grados con la función raster::terrain(). El siguiente paso será hacer coincidir espacialmente ambos raster con la misma resolución (resample()). Finalmente lo reclasificaremos teniendo en cuenta lo siguiente:

  • 1: zonas de umbría alta (Norte).

  • 2: zonas de umbría (Noreste, Noroeste).

  • 3: zonas medias (Este, Oeste).

  • 4: zonas de solana (Sureste, Suroeste).

  • 5: zonas de solana alta (Sur).

# Cargar capa y cortar a la extensión de Yunquera
mdt <- raster(here("Cartografia/CartografiaOriginal/PNOA_MDT05_ETRS89_HU30_1051_LID.tif"),
              crs = "+proj=utm +zone=30 +ellps=GRS80 +units=m +no_defs") 

# Cálculo orientaciones
orientaciones <- terrain(mdt,
                         opt = 'aspect',
                         unit = 'degrees') |> 
  resample(R_pinsapo)

# Matriz de clasificación
mat <- matrix(c(-Inf, 22.5,1,
                22.5,67.5,2,
                67.5,112.5,3,
                112.5,157.5,4,
                157.5,202.5,5,
                202.5,247.5,4,
                247.5,292.5,3,
                292.5,337.5,2,
                337.5,Inf,1),
              ncol = 3,
              byrow = TRUE)

# Reclasificación y recorte al área de estudio
R_orientaciones <- orientaciones |> 
  reclassify(rcl = mat) |> 
  mask(yunquera)


Distancia a canales

El ráster de distancia a canales se ha calculado previamente siguiendo los pasos del guión en QGIS debido a la complicidad de realizar este análisis en R. La capa obtenida se muestra en el siguiente webmap, la cual se encuentra en el sistema de coordenadas correspondiente y con resample() nos aseguramos de que coincida espacialmente con el resto de capas.

R_dist_canales <- raster(here("Cartografia/CartografiaOriginal/Raster_dist_cauces.tif")) |> 
  resample(R_pinsapo) |> 
  mask(yunquera)


Indice de Shannon

El indice de Shannon se ha calculado en la sección 2.2.3. En este caso nos quedaría rasterizar el resultado y remuestrearlo para que coincida espacialmente con el resto de capas.

R_shannon <- SF_shannon |>
  dplyr::select(shannon) |> 
  st_rasterize() |> 
  as("Raster") |> 
  resample(R_pinsapo) |> 
  mask(yunquera)


Guardar las capas

Finalmente se exportan las capas a la carpeta CapasDefRaster.

# writeRaster(R_pinsapo,
#             filename = here("Cartografia/CapasDefRaster/R_pinsapo.tif"),
#             overwrite = T)
# 
# writeRaster(R_orientaciones,
#             filename = here("Cartografia/CapasDefRaster/R_orientaciones.tif"),
#             overwrite = T)
# 
# writeRaster(R_dist_canales,
#             filename = here("Cartografia/CapasDefRaster/R_dist_canales.tif"),
#             overwrite = T)
# 
# writeRaster(R_shannon,
#             filename = here("Cartografia/CapasDefRaster/R_shannon.tif"),
#             overwrite = T)


2.3.2 Raster virtual

El ráster virtual se ha creado en QGIS mediante la herramienta Raster -> Micelanea -> Crear raster virtual.

El resultado se muestra en la Fig. 5.

<center>Fig. 5. Creación del Raster virtual en QGIS a partir de las cuatro capas de entrada</center>

Fig. 5. Creación del Raster virtual en QGIS a partir de las cuatro capas de entrada


2.3.3 Segmentación del territorio

El último paso de esta sección consiste en realizar una segmentación del monte de Yunquera utilizando el ráster virtual creado en la sección anterior. Para ello, se utilizará la herramienta Orfeo Toolbox en la interfaz de QGIS.

Para llevar esta tarea a cabo, con Orfeo correctamente instalado, se debe ir a Caja de herramientas de Procesos -> OTB -> Segmentation -> Segmentation. Una vez ahí, se seleccionan las opciones de la Fig. 6. La opción “Minimum region size” se probó con varios valores (el resultado debe convertirse a polígono). En el webmap que se verá posteriormente se analizarán las diferencias e influencia de las diferentes capas de entrada en los resultados.

<center>Fig. 6. Parámetros utilizados en la segmentación</center>

Fig. 6. Parámetros utilizados en la segmentación


Para finalizar la segunda parte de la práctica, se cargan las tres capas creadas con Orfeo en QGIS.

segmentation1500 <- read_sf(here("Cartografia/CapasDefVectorial/segmentacion1500.shp"))

segmentation2000 <- read_sf(here("Cartografia/CapasDefVectorial/segmentacion2000.shp"))

segmentation2500 <- read_sf(here("Cartografia/CapasDefVectorial/segmentacion2500.shp"))

En el webmap siguente se puede jugar con las capas para ver el efecto de las distintas segmentaciones en cada una de las capas.


Nota: los bordes de la segmentación no coinciden totalmente con el borde de los ráster. Esto se debe a que las capas se encuentran en crs 25830, pero el webmap se representa en Web Mercator (EPSG: 3857). A modo de análisis no afecta, pero las conclusiones finales se tomarán sobre el Web Map publicado de la siguiente sección.

Los resultados se discuten en la última sección. No obstante, Los resultados se harán en base a la segmentación escogida que se discute en las siguientes líneas.

Para analizar los resultados podemos activar en primer lugar Segmentacion (2000), y a continuación Segmentacion (2500) (el orden es importante). Tras hacer esto, podemos activar las variables una a una empezando por Densidad pinsapo, Indice de Shannon, Orientaciones y Distancia a canales (m). La segmentación de 2500 píxeles como unidad mínima se tratará como segmentación inicial, y las unidades segmentadas como cantones.

Lo primero que podemos ver es que si segmentamos con unidad mínima de 2000, se crean 11 cantones a mayores. La segmentación inicial parece separar muy bien áreas con alta densidad de pinsapo de áreas con baja densidad. Al aumentar el número de cantones no se ve una clara separación de nuevas áreas con densidad de pinsapo diferentes (excepto dos cantones al sur).

Si activamos el índice de Shannon vemos que no existe una clara relación entre los cantones y la diversidad de especies. Aunque hay algunos cantones que sí tienen una buena separación por diversidad, no se ve el patrón claro que se veía en la capa anterior. Además, aumentar el número de cantones no parece crear unidades homogéneas en biodiversidad.

Con respecto a la capa de orientaciones, es la capa más irregular y por ello también hace más complicado el análisis. Hay varias zonas que se ve una delimitación más o menos clara de umbrías y solanas, especialmente al noroeste del área de estudio. El aumento del número de cantones no produce claras separaciones a mayores de solanas y umbrías.

La capa de distancia a canales sí que parece tener una influenia enorme en la generación de cantones. Además, los cantones generados a mayores con segmentación de 2000 píxeles parece que más o menos todos siguen patrones de esta capa.

Si hacemos el mismo procedimiento tomando como capa inicial Segmentación (2000) y la comparamos con Segmentación (1500), vemos que los cantones generados a mayores siguien patrones sobre las capas de orientaciones y distancia a canales.

Si consideramos una gestión del pinsapo en este monte: el pinsapo es una especie que dadas sus características ecológicas se desarrolla mejor en laderas de umbría. Para realizar una gestión más eficiente de la especie, lo óptimo sería dividir el terreno en un número mayor de cantones. En este caso, con una unidad mínima de 15 hectáreas (1500 píxeles de 100 \(m^2\)) tenemos un total de 88 cantones.


2.4 Webmap

El cuarto objetivo era crear un webmap y publicarlo. Para ver dicho mapa, hacer click aquí.

2.5 Cuestiones finales

Consideras que el resultado es consistente con la realidad que conoces del terreno? Justifica tu respuesta.

Teniendo en cuenta las variables consideradas sí que es consistente el resultado. No obstante, si superponemos la segmentación junto a una imagen aérea podemos ver que al no tener en cuenta los caminos o pistas, la segmentación ignora estos elementos del paisaje que son clave en una segmentación para una gestión más sencilla y eficiente. Por ello, a efectos prácticos pierde consistencia.

Cuál de las capas (variables) utilizada aparenta tener más influencia en el resultado final? Por qué crees que pasa esto?

Distancia a canales > densidad pinsapo > orientación > indice de Shannon.